import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from matplotlib.animation import ArtistAnimation
from IPython import display
from IPython.display import HTML
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
# %matplotlib inline
# plt.rcParams['animation.ffmpeg_path'] = '/opt/homebrew/bin/ffmpeg' # It's needed to convert animation plot to html5 video
# Lax-Wendroff two step method
def lax_w_2step(u, calw):
# Firts step
w = calw(u)
u_first = np.empty_like(u)
u_first[:,0:-1] = u[:,0:-1] - kap * (w[:, 1:] - w[:,0:-1])
u_first[:, -1] = u_first[:, 0] # BdC
# Second step
w = calw(u_first)
u_second = np.empty_like(u)
u_second[:, 1:] = 0.5 * (u[:, 1:] + u_first[:, 1:] - kap * (w[:, 1:] - w[:, :-1]))
u_second[:, 0] = u_second[:, -1] # BdC
return u_second
def calw_fluid(u):
w = np.empty_like(u)
rho = u[0].copy()
rhov = u[1].copy()
e = u[2].copy()
v = rhov / rho
p = (gamma - 1) * (e - rhov * v / 2)
w[0] = rhov
w[1]= rhov * v + p
w[2] = (e + p) * v
return w
# Lax-Wendroff two step method for 2-dimention
def lax_w_2step_2d(u, calw):
# Firts step
w1, w2 = calw(u)
u_first = np.ones_like(u)
u_first[:,:-1, :-1] = u[:,:-1, :-1]\
- kap1 * (w1[:, :-1, 1:] - w1[:,:-1, :-1])\
- kap2 * (w2[:, 1:, :-1] - w2[:,:-1, :-1])
u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
# Second step
w1, w2 = calw(u_first)
u_second = np.ones_like(u)
u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:]\
- kap1 * (w1[:, 1:, 1:] - w1[:, 1:, :-1])\
- kap2 * (w2[:, 1:, 1:] - w2[:, :-1, 1:]))
u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
return u_second
def calw_fluid_q4v2(u):
w = np.empty_like(u)
rho = u[0].copy()
rhov = u[1].copy()
e = u[2].copy()
v = rhov / rho
p = (gamma - 1) * (e - rhov * v / 2)
w[0] = rhov
w[1]= rhov * v + p
w[2] = (e + p) * v
w[3] = w[2].copy()
return w, w
def calw_fluid_2d(u):
w1 = w2 = np.empty_like(u)
rho = u[0].copy()
rho_vx = u[1].copy()
rho_vy = u[2].copy()
e = u[3].copy()
vx = rho_vx / rho
vy = rho_vy / rho
v2 = vx ** 2 + vy ** 2 # velocity squared
p = (gamma - 1) * (e - rho * v2 / 2)
w1[0] = rho_vx
w1[1] = rho_vx * vx + p
w1[2] = rho_vx * vy
w1[3] = (e + p) * vx
w2[0] = rho_vy
w2[1] = rho_vy * vx
w2[2] = rho_vy * vy + p
w2[3] = (e + p) * vy
return w1, w2
差分方程式は $$ u_i^{n+1} = u_i^n + \frac{s}{h}(u_{i+1}^n - 2u_i^n + u_{i-1}^m) \\ (n = 0, 1, 2, \dots, i = 1, 2, \dots, I -1) $$
# 代入とコピーでは意味合いが違うので要注意
imax = 20
x = np.linspace(0, 1, imax + 1)
y = x.copy()
z = x
x += 2
print(f"x = \n{x}")
print(f"y = \n{y}")
print(f"z = \n{z}")
x = [2. 2.05 2.1 2.15 2.2 2.25 2.3 2.35 2.4 2.45 2.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3. ] y = [0. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1. ] z = [2. 2.05 2.1 2.15 2.2 2.25 2.3 2.35 2.4 2.45 2.5 2.55 2.6 2.65 2.7 2.75 2.8 2.85 2.9 2.95 3. ]
imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.sin(x * np.pi)
plt.plot(x, u, label=f"u0")
# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
u0 = u.copy()
u = np.array([u[i] + kap * (u[i+1] - 2*u[i] + u[i-1]) for i in range(1,imax)])
u = np.append(0, u)
u = np.append(u, 0)
if (n+1) % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x11f77d3d0>]
ちなみに、pythonではiについては次のようにfor文を用いずに計算をすることもできる。
下の関数でステップを一回進めることができる。(よって、nについてのループは必要。)
なお、bd0, bd1はx=0,1での境界条件。
# Finite-Difference Method(名前があっているかわからないけど)
def fdm(u, bd0=0, bd1=0):
w = np.empty(len(u))
w[1:-1] = u[1:-1] + kap * (u[2:] - 2*u[1:-1] + u[:-2])
w[0] = bd0
w[-1] = bd1
return w
imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.sin(x * np.pi)
plt.plot(x, u, label=f"u0")
# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
u = fdm(u)
if (n+1) % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x12f3a0c70>]
かなりすっきりとしたコードになり、再利用もできるため関数化は結構便利。
第4章例題4.3(シルクハット型関数の初期条件)
imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.zeros(imax + 1)
u[int(0.25/dx):int(0.75/dx + 1)] = 1
plt.plot(x, u, label=f"u0")
# step 3
dt = 1e-3
kap = dt/dx**2
nmax = 70
for n in range(nmax):
u = fdm(u)
if (n+1) % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x12eb79dc0>]
imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.zeros(imax + 1)
u[int(0.25/dx):int(0.75/dx + 1)] = 1
# step 3
dt = 2e-3
kap = dt/dx**2
nmax = 3
for n in range(nmax):
plt.plot(x, u, label=f"t = {dt * n:.3f}")
u = fdm(u)
plt.plot(x, u, label=f"t = {dt * (n + 1):.3f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x12f389070>]
初期条件は $$ u(x,0) = \Lambda(x) = \begin{cases} 1 - 4|x -\frac{1}{2}| & (|x - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x -\frac{1}{2}| \le \frac{1}{2}) \end{cases} $$
imax = 20
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
# step 3
dt = 1e-3
kap = dt/dx**2
# nmax = 1000
nmax = 70
fname = "data/course1.dat"
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
if n % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
np.savetxt(f, u.reshape(1, len(u))) # u[0]を末尾に書き込む
u = fdm(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.legend()
plt.xlabel("$x$")
plt.ylabel(r"$u$")
f.close()
方程式と初期、境界条件: $$ \frac{∂u}{∂t} = - \frac{∂u}{∂x}\\ u(0,t) = u(1,t) = 0 (0<t) \\ u(x,0) = \Lambda(x) \equiv \begin{cases} 1 - 4|x -\frac{1}{2}| & (|x - \frac{1}{2}| \le \frac{1}{4}) \\ 0 & (\frac{1}{4} \le |x -\frac{1}{2}| \le \frac{1}{2}) \end{cases} $$
解析解は$u(x,t) = \Lambda(x-t) \ (t \le 0.25)$である。 ラックス・ヴェンドロフ法によって差分方程式を作ると次のようになる。
先の課程の最後に作成したプログラムの差分方程式を書き換える。
imax = 50
dx = 1/imax
dt = 1e-2
kap = dt/dx
nmax = 10
fname = "data/course2.dat"
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
# step 3
for n in range(nmax):
np.savetxt(f, u.reshape(1, len(u))) # u[0]を末尾に書き込む
if n % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
u0 = u.copy()
for i in range(imax):
u[i] = u0[i] - 0.5 * kap * (u0[i+1] - u0[i-1]) + 0.5 * kap**2 * (u0[i+1] - 2*u0[i] + u0[i-1])
u[0] = 0
u[-1] = 0
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()
コードの可読性を高め、さらに改変が容易になるようにまとまった処理をするコードを関数としてひとまとめにしておくことは、いかなるプログラムにおいても有用である。 今回のコードではラックス・ヴェンドロフ法による計算で各ステップごとに同じ処理を行っているため、これを関数としてひとまとめにしておこう。
念の為境界条件も引数として準備しておくが、境界条件を変える必要がないものは特に設定する必要はない。 もしくはグローバル変数を用いて設定しても良いだろう。
# Lax-Wendroff method
def lax_wendroff(u, bd0=0, bd1=0):
w = np.empty(len(u))
w[1:-1] = u[1:-1] - 0.5 * kap * (u[2:] - u[:-2]) + 0.5 * kap ** 2 * (u[2:] - 2 * u[1:-1] + u[:-2])
w[0] = bd0
w[-1] = bd1
return w
imax = 50
dx = 1/imax
dt = 1e-2
kap = dt/dx
nmax = 10
fname = "data/course2v2.dat"
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
# step 3
for n in range(nmax):
np.savetxt(f, u.reshape(1, len(u))) # u[0]を末尾に書き込む
if (n) % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()
先のコードで確認を忘れていたが、for文も関数も結局同じであることを確認しておく。
!diff -s ./data/course2v2.dat ./data/course2.dat
Files ./data/course2v2.dat and ./data/course2.dat are identical
imax = 50
dt = 3e-2
kap = dt/dx
nmax = 10
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x + 0.1, u, label='解析解')
# step 3
for n in range(nmax):
if (n) % 10 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
# plt.ylim(-1, 1)
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x12e6d41c0>]
右にずれていく
imax = 50
dt = 2e-2
kap = dt/dx
nmax = 10
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
# step 3
for n in range(nmax):
if (n) % 5 == 0:
plt.plot(x, u, label=f"t = {dt * n:.2f}")
u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
plt.plot([0.5,0.5])
[<matplotlib.lines.Line2D at 0x11f97bca0>]
imax = 50
dt = 1e-2
dx = 1/imax
kap = dt/dx
nmax = 10
x = np.linspace(0, 1, imax + 1)
# step 2
u = np.array([1 if np.abs(x[i] - 0.5) <= 1./6. else 0 for i in range(imax+1)])
plt.plot(x, u, label=f"t = 0")
plt.plot(x + 0.1, u, label='解析解')
# step 3
for n in range(nmax):
u = lax_wendroff(u)
plt.plot(x, u, label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
<matplotlib.legend.Legend at 0x12f3c2ac0>
波形の両端で誤差が生じていることがわかる。
いよいよ2階線形波動方程式$\frac{∂^2u}{∂t^2} = \frac{∂^2 u}{∂x^2}$を数値的に解いてゆく。数学モデルは次のように、2つの一階偏微分方程式に分けられる。
ここでベクトル $$ \boldsymbol{u} = \begin{pmatrix} u \\ v \end{pmatrix} \hspace{4em} \boldsymbol{w} = \begin{pmatrix} u \\ v \end{pmatrix} $$ を用いると次のように書き改められる:
差分方程式は次のようになる:
ここで差分方程式を成分別に書くと
先のコードを二次元用に書き換える
imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course3-1.dat"
n = 0
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0]
plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u0 = u.copy()
for i in range(1, imax):
for q in range(qmax):
u[q][i] = np.array([u0[q-1][i] - 0.5 * kap * (u0[q-1][i+1] - u0[q-1][i-1]) + 0.5 * kap**2 * (u0[q][i+1] - 2*u0[q][i] + u0[q][i-1])])
for q in range(qmax):
u[q][0] = 0
u[q][-1] = 0
if not np.all(u[0] == u[1]):
print('一つ目と二つ目で値が異なります')
plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
plt.legend()
f.close()
文字の出力がないので、各ステップにおいてu[0]とu[1]が等しいことがわかる。
!diff -s ./data/course2.dat ./data/course3-1.dat
Files ./data/course2.dat and ./data/course3-1.dat are identical
結果は課程2、3で一致している
# Lax-Wendroff method for 2-d array
def lax_w_2d(u):
v = np.empty_like(u)
for q in range(qmax):
v[q][1:-1] = u[q][1:-1] - 0.5 * kap * (u[q-1][2:] - u[q-1][:-2]) + 0.5 * kap ** 2 * (u[q][2:] - 2 * u[q][1:-1] + u[q][:-2])
v[q][0] = v[q][-1] = 0
return v
imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course3-1_v2.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2d(u)
if not np.all(u[0] == u[1]):
print('一つ目と二つ目で値が異なります')
plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()
!diff -s ./data/course3-1_v2.dat ./data/course3-1.dat
Files ./data/course3-1_v2.dat and ./data/course3-1.dat are identical
(これって上で書いたコードと同じになるけど、手順分かていないのでは?)
imax = 50
nmax = 10
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-2.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.2f})")
plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()
imax = 50
nmax = 10
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-3.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = -u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.2f})")
plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()
imax = 50
nmax = 30
qmax = 2
dt = 1e-2
kap = dt/dx
dx = 1/imax
fname = "data/course3-4.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
plt.plot(x, u[0], label=f"u[0](t = 0)")
plt.plot(x, u[1], label=f"u[1](t = 0)")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
# f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
u = lax_w_2d(u)
if (n+1) % 10 == 0:
plt.plot(x, u[0], label=f"u[0](t = {dt * (n+1):.1f})")
plt.plot(x, u[1], label=f"u[1](t = {dt * (n+1):.1f})")
# np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
# f.close()
<matplotlib.legend.Legend at 0x128a7cf70>
imax = 50
nmax = 10
qmax = 2
dt = 3e-2
kap = dt/dx
dx = 1/imax
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
for n in range(nmax):
u = lax_w_2d(u)
plt.plot(x, u[0], label=f"u[0](t = {dt * (n + 1):.2f})")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
<matplotlib.legend.Legend at 0x128c45cd0>
# Lax-Wendroff method for 2-d array
def lax_w_ref(u):
v = np.empty_like(u)
for q in range(qmax):
v[q][1:-1] = u[q][1:-1] - 0.5 * kap * (u[q-1][2:] - u[q-1][:-2]) + 0.5 * kap ** 2 * (u[q][2:] - 2 * u[q][1:-1] + u[q][:-2])
v[0][0] = v[0][-1] = 0
v[1][0] = v[1][1]
v[1][-1] = v[1][-2]
return v
imax = 50
nmax = 300
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_ref(u)
txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0], "black")
ims.append([im, txt])
if (n+1) % 30 == 0:
ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig) # avoid plotting a spare static plot
これまでは線形波動方程式をラックス・ヴェンドロフ法を用いて計算してきた。これからは非線形の計算を行うために、2段階ラックス・ヴェンドロフ法を用いたコードに書き換えていく。 2段階ラックス・ヴェンドロフ法にするには差分方程式を $$ \begin{aligned} \boldsymbol{u}_i^{\overline{n+1}} &= \boldsymbol{u}_i^n - \frac{s}{h}(\boldsymbol{w}_{i+1}^n - \boldsymbol{w}_i^n) \\ \boldsymbol{u}_i^{n+1} &= \frac{1}{2}\left[\boldsymbol{u}_i^n + \boldsymbol{u}_i^{\overline{n+1}} - \frac{s}{h}(\boldsymbol{w}_{i}^{\overline{n+1}} - \boldsymbol{w}_{i-1}^{\overline{n+1}})\right] \end{aligned} $$ とすればよい。 ここで新たな変数$\boldsymbol{w}_1^n$と$\boldsymbol{u}_i^{\overline{n+1}}$が必要となる。
imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-1.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappendモードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u0 = u.copy()
w[1] = u[0].copy()
w[0] = u[1].copy()
for q in range(qmax):
u[q][0] = 0
u[q][-1] = 0
for i in range(1, imax):
u[q][i] = np.array([u0[q][i] - 0.5 * kap * (w[q][i+1] - w[q][i-1]) + 0.5 * kap**2 * (u0[q][i+1] - 2*u0[q][i] + u0[q][i-1])])
plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()
ここでも結果に変化がないことを確認しておく
! diff -s ./data/course2.dat ./data/course4-1.dat
Files ./data/course2.dat and ./data/course4-1.dat are identical
imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-2.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u0 = u.copy()
w[1] = u[0].copy()
w[0] = u[1].copy()
for i in range(0, imax):
for q in range(qmax):
uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
uh[0][imax] = 0.0
uh[1][imax] = uh[1][imax - 1]
w[0] = uh[1].copy()
w[1] = uh[0].copy()
for i in range(1, imax + 1):
for q in range(qmax):
u[q][i] = 0.5 * (u0[q][i] + uh[q][i] - kap * (w[q][i] - w[q][i-1]))
u[0][0] = 0.0
u[1][0] = u[1][1]
plt.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u$")
plt.legend()
f.close()
今回も積分の計算は関数としてひとまとめにしておこう。
関数化するにあたって、wの計算が必要になるため、それをcalw()として引数にとることにする。
また、先の関数では境界条件を引数に設定していたが、ここでは手順4で課す周期的境界条件を課しておくことにする。
# Lax-Wendroff two step method
def lax_w_2step(u, calw):
# Firts step
w = calw(u)
u_first = np.empty_like(u)
u_first[:,0:-1] = u[:,0:-1] - kap * (w[:, 1:] - w[:,0:-1])
u_first[:, -1] = u_first[:, 0] # BdC
# Second step
w = calw(u_first)
u_second = np.empty_like(u)
u_second[:, 1:] = 0.5 * (u[:, 1:] + u_first[:, 1:] - kap * (w[:, 1:] - w[:, :-1]))
u_second[:, 0] = u_second[:, -1] # BdC
return u_second
実習資料に合わせるならば
def calw(u, w):
w[0] = u[1].copy()
w[1] = u[0].copy()
とすれば良いだろう。
しかし、このような関数内部で変数を変更しようとすると、下の「おまけ」で見るように代入するもの(配列、文字列、配列の要素、etc.)によって処理が変わってしまうので、仕様を理解していないと意図しない挙動になりそうである。 (配列の要素の代入はきちんと置き換えられるため上の関数でも良いのだが)
そこで、次のように0、1要素を入れ替えた配列(ndarray)を返す関数を作り、wに代入するほうが良いのではないだろうか。
def calw_wave(u):
w = np.empty_like(u)
w[0] = u[1].copy()
w[1] = u[0].copy()
return w
# example of usage
# w = calw_wave(u)
この関数の実行例を見て見ると、きちんと第0要素と第1要素が入れ替わっていることがわかる。
# Test
u = np.arange(0,10).reshape(2, 5)
w = calw_wave(u)
print(f"u = \n{u}")
print(f"w = \n{w}")
print(f"u[0] == w[1]: {np.all(u[0] == w[1])}")
print(f"u[1] == w[0]: {np.all(u[1] == w[0])}")
u = [[0 1 2 3 4] [5 6 7 8 9]] w = [[5 6 7 8 9] [0 1 2 3 4]] u[0] == w[1]: True u[1] == w[0]: True
ndarrayの要素に代入(変更される)
def f(x,y):
for i in range(10):
y[i] = x[i]
print(f"関数内:{y}")
x = np.arange(0,10,1)
y = np.zeros(10)
print(f"初め:{y}")
f(x, y)
print(f"関数実行後:{y}")
初め:[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] 関数内:[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] 関数実行後:[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
ndarray全体を代入(変更されない)
def f(x,y):
y = x.copy()
print(f"関数内:{y}")
x = np.arange(0,10,1)
y = np.zeros(10)
print(f"初め:{y}")
f(x, y)
print(f"関数実行後:{y}")
初め:[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] 関数内:[0 1 2 3 4 5 6 7 8 9] 関数実行後:[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
文字列(変更されない)
def f(a):
a += "def"
print(f"関数内:{a}")
a = "abc"
print(f"初め:{a}")
f(a)
print(f"関数実行後:{a}")
初め:abc 関数内:abcdef 関数実行後:abc
多次元のndarrayの要素の配列を代入(変更される)
def f(x, y):
y[0] = x[1].copy()
y[1] = x[0].copy()
x = np.arange(0,10).reshape(2,5)
y = np.zeros(10).reshape(2,5)
print(f"初め: y =\n{y}")
f(x, y)
print(f"関数実行後: y = \n{y}")
初め: y = [[0. 0. 0. 0. 0.] [0. 0. 0. 0. 0.]] 関数実行後: y = [[5. 6. 7. 8. 9.] [0. 1. 2. 3. 4.]]
imax = 50
nmax = 10
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-2_v2.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f" t = 0")
plt.plot(x + 0.1, u[0], label='解析解')
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_wave)
plt.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
f.close()
!diff -s ./data/course4-2_v2.dat ./data/course4-2.dat
Files ./data/course4-2_v2.dat and ./data/course4-2.dat are identical
今の計算範囲では境界条件が使われていないため、関数版(*_v2.dat)とfor文版で結果が同じであることがわかる。
imax = 50
nmax = 50
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
plt.plot(x, u[0], label=f"t = 0")
for n in range(nmax):
u0 = u.copy()
w[1] = u[0].copy()
w[0] = u[1].copy()
for i in range(0, imax):
for q in range(qmax):
uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
uh[0][imax] = 0.0
uh[1][imax] = uh[1][imax - 1]
w[0] = uh[1].copy()
w[1] = uh[0].copy()
for i in range(1, imax + 1):
for q in range(qmax):
u[q][i] = 0.5 * (u0[q][i] + uh[q][i]) - kap * (w[q][i] - w[q][i-1])
u[0][0] = 0.0
u[1][0] = u[1][1]
plt.plot(x, u[0], label=f"t = {dt * (n + 1):.2f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend()
<matplotlib.legend.Legend at 0x12f971880>
関数の書き換えの手間を省くため、ここでは関数化前のコードを利用した。
形が崩れていることがわかる。
一周するようにnmaxを100に変えておく。
imax = 50
nmax = 100
qmax = 2
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course4-4.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
# plt.plot(x, u[0], label=f"t = 0")
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_wave)
txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0], "black")
ims.append([im, txt])
if (n+1) % 30 == 0:
ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig) # avoid plotting a spare static plot
いよいよ非線形方程式の数値解法に入る。
qmaxを3に増やす¶まず、qmaxの値を3に増やし、増やした2番目の要素は1番目の要素と同じになるように書き換える。
つまり
u[2][i] = u[1][i]
となるようにする。
この時、先に定義したcalw_wave()では対応できないため、これも少し修正する。(もっとも、ここでしか使わないため適当である。)
def calw_wave_q3(u):
w = np.empty_like(u)
w[0] = u[1].copy()
w[1] = u[0].copy()
w[2] = w[1].copy()
return w
imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
fname = "data/course5-1.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.array([1 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0 for i in range(imax + 1)])
u[1] = u[0].copy()
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
if (n) % 20 == 0:
plt.plot(x, u[0], label=f"t = {dt * n:.2f}")
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_wave_q3)
plt.plot(x, u[0], label=f"t = {dt * (n+1):.1f}")
plt.xlabel("$x$")
plt.ylabel(r"$u[0]$")
plt.legend(loc="upper left")
f.close()
!diff -s ./data/course4-4.dat ./data/course5-1.dat
Files ./data/course4-4.dat and ./data/course5-1.dat are identical
計算結果が先のものと変わらないことも確かめられた。
流速密度行列$\boldsymbol{w}$の計算を $$ \boldsymbol{w} = \begin{pmatrix} \rho v \\ \rho v^2 + p \\ (e + p) v \end{pmatrix} $$ に対応するように書き換える。ここで、$\rho$は気体の密度、$v$は気体の速度、$e$はエネルギー密度である。 この時、各メッシュにおいて $$ \boldsymbol{u} = \begin{pmatrix} \rho \\ \rho v \\ e \end{pmatrix} $$ の成分を用いて $$ u = \frac{\rho v}{\rho}, \hspace{4em} p = (\gamma - 1)\left(e - \frac{\rho}{2}v^2 \right) $$ と計算する。
def calw_fluid(u):
w = np.empty_like(u)
rho = u[0].copy()
rhov = u[1].copy()
e = u[2].copy()
v = rhov / rho
p = (gamma - 1) * (e - rhov * v / 2)
w[0] = rhov
w[1]= rhov * v + p
w[2] = (e + p) * v
return w
imax = 50
nmax = 50
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.ones(imax+1)
u[1] = np.zeros(imax+1)
u[2] = np.array([1.9 - 4 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0.9 for i in range(imax + 1)])
fig, axes = plt.subplots(1,3,figsize=(12, 3)) # for statical plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4)) # for animation plot
ims = []
plots = [0] * 3
for i, ax in enumerate(axes):
ax.plot(x, u[i], label="t = 0")
for n in range(nmax):
# animation
for i, ax_anim in enumerate(axes_anim):
plots[i] = ax_anim.plot(x, u[i], c="black")
txt = axes_anim[0].text(0, 1.5, f"t = {dt * (n):.1f}", size=16)
plots[0] = [plots[0][0],txt] # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2]) # for subplots in animation
u = lax_w_2step(u, calw_fluid)
# nomal plot
if (n + 1) % 10 == 0:
for i, ax in enumerate(axes):
ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")
# Set x, and y labels
for axs in [axes, axes_anim]:
axs[1].set_xlabel("$x$")
axs[0].set_ylabel(r"$\rho$")
axs[1].set_ylabel(r"$\rho v$")
axs[2].set_ylabel(r"energy")
axes[0].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=100)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim) # close animation plot
A = 0.001¶imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 1e-3
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = 1 + A * np.cos(2*np.pi*x)
u[1] = A * np.cos(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.cos(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[0]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_fluid) # 積分1ステップ進める
# animation
txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0], "black")
ims.append([im, txt])
if (n+1) % 25 == 0:
ax2.plot(x, u[0], ".", label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
この時には安定した波動が観察される。
A = 0.1¶この時は非線形の影響が強く出ることになる。
imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/course5-ex1-A01.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[1]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0].reshape(1, len(u[1]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_fluid) # 積分1ステップ進める
# animation
txt = ax.text(0, 0.3, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0], "black")
ims.append([im, txt])
if (n+1) % 100 == 0:
ax2.plot(x, u[0], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.
方程式の非線形な部分が効いており、衝撃波のようなものが観察される。
初期設定が、実習書では$\cos$、教科書では$\sin$が使われていたが、それに気づかず実行結果が異なって見えたため、関数化の際の間違いがないかを確認するべくfor文を回して確認した。
実際は初期設定を変えることで解決したが、関数化の妥当性を見るために以下の記述を残しておく。
# for文を使った計算(検算用)
imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/course5-ex1-A01_check.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
# plt.plot(x, u[0], label=f"t = 0")
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[1], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[1]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[1].reshape(1, len(u[1]))) # u[1]を末尾に書き込む
u0 = u.copy()
w = calw_fluid(u)
for i in range(0, imax):
for q in range(qmax):
uh[q][i] = u0[q][i] - kap * (w[q][i+1] - w[q][i])
uh[:, -1] = uh[:,0]
w = calw_fluid(uh)
for i in range(1, imax + 1):
for q in range(qmax):
u[q][i] = 0.5 * (u0[q][i] + uh[q][i] - kap * (w[q][i] - w[q][i-1]))
u[:, 0] = u[:, -1]
# animation
txt = ax.text(0, 0.95, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[1], "black")
ims.append([im, txt])
if (n+1) % 100 == 0:
ax2.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.
!diff -s ./data/course5-ex1-A01_check.dat ./data/course5-ex1-A01.dat
Files ./data/course5-ex1-A01_check.dat and ./data/course5-ex1-A01.dat are identical
imax = 50
nmax = 100
qmax = 3
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = np.ones(imax+1)
u[1] = np.zeros(imax+1)
u[2] = np.array([0.4 + 2 * np.abs(x[i] - 0.5) if np.abs(x[i] - 0.5) <= 0.25 else 0.9 for i in range(imax + 1)])
fig, axes = plt.subplots(1,3,figsize=(12, 3)) # for static plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4)) # for animation plot
ims = []
plots = [0] * 3
labels = [r"$\rho$", r"$\rho v$", r"energy"]
# set first plot and labesl
for i, (ax, ax_anim) in enumerate(zip(axes, axes_anim)):
ax.plot(x, u[i], label="t = 0")
plots[i] = ax_anim.plot(x, u[i], c="black")
ax.set_ylabel(labels[i])
ax_anim.set_ylabel(labels[i])
ax.set_xlabel(r"$x$")
ax_anim.set_xlabel(r"$x$")
txt = axes_anim[0].text(0, 1.8, f"t = 0", size=16)
plots[0] = [plots[0][0],txt] # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2]) # for subplots in animation
for n in range(nmax):
u = lax_w_2step(u, calw_fluid)
# animation
for i, ax_anim in enumerate(axes_anim):
plots[i] = ax_anim.plot(x, u[i], c="black")
txt = axes_anim[0].text(0, 1.5, f"t = {dt * (n):.1f}", size=16)
plots[0] = [plots[0][0],txt] # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2]) # for subplots in animation
# nomal plot
if (n + 1) % 50 == 0:
for i, ax in enumerate(axes):
ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")
axes[1].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=100)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim) # close animation plot
imax = 50
nmax = 200
qmax = 3
dx = 1/imax
dt = 1e-3
kap = dt/dx
gamma = 5/3
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
rho = np.ones(imax+1)
v = np.zeros(imax+1)
p = 0.9 * np.ones_like(x)
rho[round(1/6/dx):round(1/3/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)]))
rho[round(2/3/dx):round(5/6/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)]))
v[round(1/6/dx):round(1/3/dx)] = -np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)])
v[round(2/3/dx):round(5/6/dx)] = -np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)])
u[0] = rho
u[1] = rho * v
u[2] = p/(gamma-1) + rho * v * v / 2
fig, axes = plt.subplots(1,3,figsize=(12, 3)) # for static plot
fig_anim, axes_anim = plt.subplots(1,3,figsize=(15, 4)) # for animation plot
ims = []
plots = [0] * 3
labels = [r"$\rho$", r"$\rho v$", r"energy"]
# set first plot and labesl
for i, (ax, ax_anim) in enumerate(zip(axes, axes_anim)):
ax.plot(x, u[i], label="t = 0")
plots[i] = ax_anim.plot(x, u[i], c="black")
ax.set_ylabel(labels[i])
ax_anim.set_ylabel(labels[i])
ax.set_xlabel(r"$x$")
ax_anim.set_xlabel(r"$x$")
txt = axes_anim[0].text(0, 1.8, f"t = 0", size=16)
plots[0] = [plots[0][0],txt] # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2]) # for subplots in animation
for n in range(nmax):
u = lax_w_2step(u, calw_fluid)
# animation
for i, ax_anim in enumerate(axes_anim):
plots[i] = ax_anim.plot(x, u[i], c="black")
txt = axes_anim[0].text(0, 1.8, f"t = {dt * (n):.1f}", size=16)
plots[0] = [plots[0][0],txt] # to add text in the first ax, this is needed
ims.append(plots[0]+plots[1]+plots[2]) # for subplots in animation
# nomal plot
if (n + 1) % 100 == 0:
for i, ax in enumerate(axes):
ax.plot(x, u[i], label=f"t = {dt * (n+1):.1f}")
axes[1].legend(loc="best")
# Adjust layouts
fig.tight_layout()
fig_anim.tight_layout()
# Start animation
anim = ArtistAnimation(fig_anim, ims, interval=10)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
plt.close(fig_anim) # close animation plot
結果がおかしいので、wの計算が積分がおかしいかどうかfor文を用いたものと比較してみる。
したのをやってみて思ったが、積分の計算は同じ結果になるのは当たり前だと思う。 wの計算も間違えようがないとも思うし、一体いどうして計算が破綻するのだろうか。
# Lax-Wendroff two step method for ver.
def lax_w_2step_for(u, calw):
# Firts step
w = calw(u)
u_first = np.empty_like(u)
for i in range(imax):
u_first[:,i] = u[:,i] - kap * (w[:,i+1] - w[:, i])
u_first[:, -1] = u_first[:, 0] # BdC
# Second step
w = calw(u_first)
u_second = np.empty_like(u)
for i in range(1, imax+1):
u_second[:, i] = 0.5 * (u[:, i] + u_first[:, i] - kap * (w[:, i] - w[:, i-1]))
u_second[:, 0] = u_second[:, -1] # BdC
return u_second
imax = 50
nmax = 350
qmax = 3
dx = 1/imax
dt = 1e-4
kap = dt/dx
gamma = 5/3
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
rho = np.ones(imax+1)
v = np.ones(imax+1)
p = np.array([0.9]*(imax+1))
rho[round(1/6/dx):round(1/3/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)]))
rho[round(2/3/dx):round(5/6/dx)] = 1 + np.abs(np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)]))
v[round(1/6/dx):round(1/3/dx)] = -np.sin(6*np.pi*x[round(1/6/dx):round(1/3/dx)])
v[round(2/3/dx):round(5/6/dx)] = -np.sin(6*np.pi*x[round(2/3/dx):round(5/6/dx)])
u[0] = rho
u[1] = rho * v
u[2] = p/(gamma-1) + rho * v * v / 2
v = u.copy()
# plt.plot(x, u[2], label=f"t = 0")
for n in range(nmax):
u = lax_w_2step(u, calw_fluid)
v = lax_w_2step_for(v, calw_fluid)
if np.any(v != u):
print("u is not equal tot\ v!!!")
今までは1次元の問題のみを扱ってきたが、これからは多次元の問題に取り組んでいく。初めは2次元である。
2次元流体力学方程式は次のように表せる。
はエネルギー密度
def calw_fluid_q4(u):
w = np.empty_like(u)
rho = u[0].copy()
rhov = u[1].copy()
e = u[2].copy()
v = rhov / rho
p = (gamma - 1) * (e - rhov * v / 2)
w[0] = rhov
w[1]= rhov * v + p
w[2] = (e + p) * v
w[3] = w[2].copy()
return w
imax = 50
nmax = 200
qmax = 4
dx = 1/imax
dt = 1e-2
kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-1-A03.dat"
x = np.linspace(0, 1, imax + 1)
u = np.zeros((qmax, imax+1))
w = np.zeros((qmax, imax+1))
uh = np.zeros((qmax, imax+1))
u[0] = 1 + A * np.sin(2*np.pi*x)
u[1] = A * np.sin(2*np.pi*x)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*x))
u[3] = u[2].copy()
# plt.plot(x, u[0], label=f"t = 0")
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[1], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[1]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[1].reshape(1, len(u[1]))) # u[0]を末尾に書き込む
u = lax_w_2step(u, calw_fluid_q4) # 積分1ステップ進める
# animation
txt = ax.text(0, 0.3, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[1], "black")
ims.append([im, txt])
if (n+1) % 100 == 0:
ax2.plot(x, u[1], label=f"t = {dt * (n+1):.2f}")
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
Font 'rm' does not have a glyph for '\uff11' [U+ff11], substituting with a dummy symbol.
!diff -s data/course5-ex1-A03.dat data/courseB-1-A03.dat
Files data/course5-ex1-A03.dat and data/courseB-1-A03.dat are identical
今回の書き換えではu[0]に変化がない。
新たな次元が加わることで、ラックス・ウェンドロフ法の関数も書き換えなければならない。 変更は以下の要領で行う。
(第1段階,予測子)
uh[q][i][j] = uo[q][i][j] - kap1*(w1[q][i+1][j]-w1[q][i][j]) - kap2*(w2[q][i][j+1]-w2[q][i][j])
(第2段階,修正子)
u[q][i][j] = 0.5*(uo[q][i][j] + uh[q][i][j] - kap1*(w1[q][i][j]-w1[q][i-1][j])
- kap2*(w2[q][i][j]-w2[q][i][j-1]))
# Lax-Wendroff two step method for 2-dimention
def lax_w_2step_2d(u, calw):
# Firts step
w1, w2 = calw(u)
u_first = np.ones_like(u)
u_first[:,:-1, :-1] = u[:,:-1, :-1] - kap1 * (w1[:, 1:, :-1] - w1[:,:-1, :-1]) - kap2 * (w2[:, :-1, 1:] - w2[:,:-1, :-1])
u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
# Second step
w1, w2 = calw(u_first)
u_second = np.ones_like(u)
u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:] - kap1 * (w1[:, 1:, 1:] - w1[:, :-1, 1:]) - kap2 * (w2[:, 1:, 1:] - w2[:, 1:, :-1]))
u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
return u_second
yに対応するu,wの次元を増やす。¶def calw_fluid_q4v2(u):
w = np.empty_like(u)
rho = u[0].copy()
rhov = u[1].copy()
e = u[2].copy()
v = rhov / rho
p = (gamma - 1) * (e - rhov * v / 2)
w[0] = rhov
w[1]= rhov * v + p
w[2] = (e + p) * v
w[3] = w[2].copy()
return w, w
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
# kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-2-A01.dat"
j_flag = 0
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1) # y軸方向の
u = np.zeros((qmax, jmax+1, imax+1)) # u[q,y,x] 次元を増やすため、jmax+1 の部分を足した。
X, Y = np.meshgrid(x, y)
u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel(r"$u[0]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[0]$")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[1, 0, :].reshape(1, len(u[1, 0, :])))
u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
# animation
txt = ax.text(0, 1, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0,0, :], "black")
ims.append([im, txt])
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
!diff -s data/courseB-2-A01.dat data/course5-ex1-A01.dat
from mpl_toolkits.mplot3d import Axes3D
としてmatplotlib の3Dプロットを用いる。
しかし、描画にものすごく時間がかかるので、3Dプロットは今後使わず、等高線プロットにしておこう。
imax = 50
jmax = 50
nmax = 50
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1) # y軸方向の
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()
# fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
fig3, ax3 = plt.subplots(figsize=(9,6),subplot_kw=dict(projection='3d')) # for 3d animation
ax2.plot(x, u[1,0,:], label="t = 0")
ims = []
# ax.set_xlabel("$x$")
# ax.set_ylabel(r"$u[1]$")
ax2.set_xlabel("$x$")
ax2.set_ylabel(r"$u[1]$")
ax3.set_xlabel("$x$")
ax3.set_ylabel("$y$")
ax3.set_zlabel("u[1]")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[1, 0, :], label=f"t = {dt * (n+1):.2f}")
txt = ax3.text(0, 0, 0, f"t = {dt * (n + 1):.1f}")
im = ax3.plot_surface(X, Y, u[0, :, :], cmap=cm.viridis)
ims.append([im, txt])
anim = ArtistAnimation(fig3, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig3)
等高線プロットはax.contourで作成する。
これはこれで時間がかかるけれど、3次元プロットよりも幾分かマシである。
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
u[3] = u[2].copy()
fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_q4v2) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video())
ax2.legend()
plt.close(fig)
calwを変える¶はエネルギー密度
def calw_fluid_2d(u):
w1 = w2 = np.empty_like(u)
rho = u[0].copy()
rho_vx = u[1].copy()
rho_vy = u[2].copy()
e = u[3].copy()
vx = rho_vx / rho
vy = rho_vy / rho
v2 = vx ** 2 + vy ** 2 # velocity squared
p = (gamma - 1) * (e - rho * v2 / 2)
w1[0] = rho_vx
w1[1] = rho_vx * vx + p
w1[2] = rho_vx * vy
w1[3] = (e + p) * vx
w2[0] = rho_vy
w2[1] = rho_vy * vx
w2[2] = rho_vy * vy + p
w2[3] = (e + p) * vy
return w1, w2
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
# kap = dt/dx
gamma = 5/3
A = 0.1
fname = "data/courseB-6-A01.dat"
j_flag = 0
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
u = np.zeros((qmax, jmax+1, imax+1)) # u[q,y,x] 次元を増やすため、jmax+1 の部分を足した。
X, Y = np.meshgrid(x, y)
u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
fig, ax = plt.subplots(figsize=(6,4))
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0, 0,:], label="t = 0")
ims = []
ax.set_xlabel("$x$")
ax.set_ylabel("u[0]")
ax2.set_xlabel("$x$")
ax2.set_ylabel("u[0]")
# テキストファイルが存在していたら削除する
if os.path.exists(fname):
os.remove(fname)
f = open(fname, "a") # 書き込み用のファイルをappend(追記)モードで開く
for n in range(nmax):
np.savetxt(f, u[0, 0, :].reshape(1, len(u[1, 0, :])))
u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
# animation
txt = ax.text(0, 1, f"t = {dt * (n + 1):.1f}")
im, = ax.plot(x, u[0, 0, :], "black")
ims.append([im, txt])
f.close()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
ax2.legend()
plt.close(fig)
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
u[0] = 1 + A * np.sin(2*np.pi*X)
u[1] = A * np.sin(2*np.pi*X)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X))
fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0, 0 ,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$y$")
ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, 0 ,:], label=f"t = {dt * (n+1):.2f}")
txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
im = ax.contourf(X, Y, u[0, : ,:], 20, cmap=plt.cm.viridis)
ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video())
ax2.legend()
plt.close(fig)
!diff -s data/courseB-6-A01.dat data/courseB-2-A01.dat
課程5のものとプロットを比較してみよう。 概ね等しいが、僅かに異なっていることがわかる。
f1 = np.loadtxt("data/course5-ex1-A01.dat")
f2 = np.loadtxt("data/courseB-6-A01.dat")
fig, ax = plt.subplots()
ax.plot(x, f1[0,:] - f2[0,:], label="n = 0")
# ax.plot(x, f1[0,:], "--", label="課程5")
# ax.plot(x, f2[0,:], ".", label="付録B")
for i in range(99, 202, 100):
# ax.plot(x, f1[i,:], "--", label="課程5")
# ax.plot(x, f2[i,:], ".", label="付録B")
ax.plot(x, f1[i,:] - f2[i,:], label=f"n = {i}")
ax.legend()
<matplotlib.legend.Legend at 0x13fe164f0>
u[i]が$\text{imax} \times \text{jmax}$行列で、$\mathbf{x} \otimes \mathbf{y}$に等しいので、それの転置を取れば$\mathbf{y} \otimes \mathbf{x}$になり、$x, y$を入れ替えた場合を考えられる。
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
u[0] = 1 + A * np.sin(2*np.pi*Y)
u[1] = A * np.sin(2*np.pi*Y)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*Y))
fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,:,0], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$y$")
ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, :, 0], label=f"t = {dt * (n+1):.2f}")
txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=30)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video())
ax2.legend()
plt.close(fig)
animationから明らかなように、y方向への伝播が観測できる。
次のコードを実行してわかるように、$t = 0.6$程度で計算が破綻している。 x,y方向への伝播ではみえなかったが、後の爆発のシムレーションの際に左上、右下の領域での計算結果がおかしくなっていることが見えるため、この非対称性によって計算がうまくいっていないのだろう。
計算が対称的でないのはラックスウェンドロフ法の関数に起因するであろうから、見直すべきかもしれない。
imax = 50
jmax = 50
nmax = 70
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-2
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
X_roll = np.zeros_like(X)
for j in range(jmax):
X_roll[j, :] = np.roll(X[j], -j)
u[0] = 1 + A * np.sin(2*np.pi*X_roll)
u[1] = A * np.sin(2*np.pi*X_roll)
u[2] = u[1].copy()
u[3] = 0.9 * (1 + gamma * A * np.sin(2*np.pi*X_roll))
# u = np.roll(u, 35, axis=1)
fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(figsize=(6,4))
ax2.plot(x, u[0,0,:], label="t = 0")
ims = []
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
if (n+1) % 100 == 0:
ax2.plot(x, u[0, 0, :], label=f"t = {dt * (n+1):.2f}")
txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
cbar = fig.colorbar(im)
cbar.set_label(r"Density $\rho$", fontsize=15)
fig.tight_layout()
anim = ArtistAnimation(fig, ims, interval=50)
video = anim.to_html5_video()
html = display.HTML(video)
display.display(html)
# HTML(anim.to_html5_video())
ax2.legend()
plt.close(fig)
X_roll = np.zeros_like(X)
for j in range(jmax):
X_roll[j, :] = np.roll(X[j], j)
plt.contourf(X, Y, X_roll)
# X[0, :]
<matplotlib.contour.QuadContourSet at 0x1681254c0>
以下、摂動を加える必要あり
レイリーテイラー不安定性 密度を上部、下部で変えて重力を設定。
ケルビンヘルムホルツ不安定性 vxを、下半分は右側、上半分は左側に動くように設定する。
imax = 50
jmax = 50
nmax = 160
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 5e-3
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
nplot = 5
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
w0 = 0.1
sigma = 0.05/np.sqrt(2.)
gamma = 5/3.
rho = np.ones_like(X)
vx = -0.5 + (np.abs(Y-0.5)<0.25)
vy = w0*np.sin(4*np.pi*X) * ( np.exp(-(Y-0.25)**2/(2 * sigma**2)) + np.exp(-(Y-0.75)**2/(2*sigma**2)) )
vz = 0*X
P = 0*X + 2.5
u[0] = rho
u[1] = rho * vx
u[2] = rho * vy
u[3] = rho * (vx ** 2 + vy ** 2) / 2 + P/(gamma - 1)
# plt.contourf(X, Y, u[1])
# plt.colorbar()
# plt.show()
# plt.contourf(X, Y, u[0])
# plt.colorbar()
# plt.show()
# u[0] = np.tile(1 + A * np.sin(2*np.pi*x), (jmax+1,1)) # density rho
# u[1] = np.tile(A * np.sin(2*np.pi*x), (jmax+1,1)) # momentum rho v
# u[2] = np.tile(0.9 * (1 + gamma * A * np.sin(2*np.pi*x)), (jmax+1,1)) # energy e
# u[3] = u[2].copy()
# fig, ax = plt.subplots(figsize=(6,4), dpi=100)
fig2, ax2 = plt.subplots(2, 1, figsize=(6,12))
# fig2.colorbar(plot)
# ims = []
# ax.set_xlabel(r"$x$")
# ax.set_ylabel(r"$y$")
# ax2.set_xlabel(r"$y$")
# ax2.set_ylabel(r"$u[0]$")
plot = ax2[0].contourf(X, Y, u[1, :, ])
# txt = ax.text(0, 1.01, f"t = 0", fontsize=15)
# im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
# ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
# if n == 0:
# ax2.contourf(X, Y, u[1, :, :])
# txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
# im = ax.contourf(X, Y, u[0, :, :], 20, cmap=plt.cm.viridis)
# ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
plot = ax2[1].contourf(X, Y, u[1, :, ])
# cbar = fig.colorbar(im)
# cbar.set_label(r"Density $\rho$", fontsize=15)
# fig.tight_layout()
# anim = ArtistAnimation(fig, ims, interval=30)
# video = anim.to_html5_video()
# html = display.HTML(video)
# display.display(html)
# ax2.legend()
# plt.close(fig)
imax = 50
jmax = 50
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 2e-3
t_end = 2
nmax = int(t_end / dt)
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
plot_index = 1
nplot = 10
mabiki = 1
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
def lax_w_2step_2d_ver2(u, calw):# 境界条件を変える
# Firts step
w1, w2 = calw(u)
u_first = np.ones_like(u)
u_first[:,:-1, :-1] = u[:,:-1, :-1] - kap1 * (w1[:, 1:, :-1] - w1[:,:-1, :-1]) - kap2 * (w2[:, :-1, 1:] - w2[:,:-1, :-1])
u_first[:, :, -1] = u_first[:, :, 0].copy() # BdC
u_first[:, -1, :] = u_first[:, 0, :].copy() # BdC
# Second step
w1, w2 = calw(u_first)
u_second = np.ones_like(u)
u_second[:, 1:, 1:] = 0.5 * (u[:, 1:, 1:] + u_first[:, 1:, 1:] - kap1 * (w1[:, 1:, 1:] - w1[:, :-1, 1:]) - kap2 * (w2[:, 1:, 1:] - w2[:, 1:, :-1]))
u_second[:, :, 0] = u_second[:, :, -1].copy() # BdC
u_second[:, 0, :] = u_second[:, -1, :].copy() # BdC
return u_second
w0 = 0.1
sigma = 0.05/np.sqrt(2.)
gamma = 5/3
delta = 0.05
rho = np.ones_like(X)
vx = np.tanh(- (np.abs(Y - 0.5) - 0.25) / delta)
vy = w0*np.sin(2*np.pi*X) * vx * np.tanh(Y / delta)
P = 0*X + 2.5
# w0 = 0.1
# sigma = 0.05/np.sqrt(2.)
# rho = 1. + (np.abs(Y-0.5) < 0.25)
# vx = -0.5 + (np.abs(Y-0.5)<0.25)
# vy = w0*np.sin(4*np.pi*X) * ( np.exp(-(Y-0.25)**2/(2 * sigma**2)) + np.exp(-(Y-0.75)**2/(2*sigma**2)) )
# P = 2.5 * np.ones(X.shape)
u[0] = rho
u[1] = rho * vx
u[2] = rho * vy
u[3] = rho * (vx ** 2 + vy ** 2) / 2 + P/(gamma - 1)
for ui in u:
plt.contourf(X, Y, ui)
plt.colorbar()
plt.show()
def plot_vector():
# plt.streamplot(X[::mabiki, ::mabiki], Y[::mabiki, ::mabiki], vx[::mabiki, ::mabiki], vy[::mabiki, ::mabiki], density=1.5)
plt.quiver(X[::mabiki, ::mabiki], Y[::mabiki, ::mabiki], vx[::mabiki, ::mabiki], vy[::mabiki, ::mabiki])
plt.title(f"n = {n}, t = {n * dt}")
plt.show()
for n in range(nmax):
if n % (nmax / nplot) == 0:
plot_vector()
u = lax_w_2step_2d_ver2(u, calw_fluid_2d) # 積分1ステップ進める
vx = u[1] / u[0]
vy = u[2] / u[0]
plot_vector()
vx = np.tanh(- (np.abs(Y - 0.5) - 0.01) / delta*10)
plt.contourf(X, Y, vx)
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x13bd32af0>
imax = 50
jmax = 50
nmax = 200
qmax = 4
dx = 1/imax
dy = 1/jmax
dt = 1e-3
kap1 = dt/dx
kap2 = dt/dy
gamma = 5/3
A = 0.1
nindex = 3
nplot = 10
x = np.linspace(0, 1, imax + 1)
y = np.linspace(0, 1, jmax + 1)
X, Y = np.meshgrid(x, y)
u = np.zeros((qmax, imax+1, jmax+1)) # 次元を増やすため、jmax+1 の部分を足した。
R = np.sqrt((X - 0.5) ** 2 + (Y - 0.5) ** 2)
u[0] = np.ones_like(X)
u[1] = np.zeros_like(X)
u[2] = np.zeros_like(X)
u[3] = 0.9 * np.ones_like(X) + (1 - 2 * R) * (R < 0.5)
for ui in u:
plt.contourf(X, Y, ui)
plt.colorbar()
plt.show()
# fig, ax = plt.subplots(figsize=(5,4), dpi=100)
# fig2, ax2 = plt.subplots(figsize=(6,4))
# ax2.plot(x, u[nindex,0,:], label="t = 0")
# ims = []
# ax.set_xlabel(r"$x$")
# ax.set_ylabel(r"$y$")
# ax.set_aspect('equal', adjustable='box')
# ax2.set_xlabel(r"$x$")
# ax2.set_ylabel(r"$u[0]$")
for n in range(nmax):
u = lax_w_2step_2d(u, calw_fluid_2d)
if (n + 1) % (nmax / nplot) == 0:
fig, ax = plt.subplots(figsize=(5,4))
ax.set_aspect('equal', adjustable='box')
ax.contourf(X, Y, u[3])
fig.suptitle(f"t = {dt * (n+1):.3f}")
fig.tight_layout()
fig.show()
# for n in range(nmax):
# u = lax_w_2step_2d(u, calw_fluid_2d) # 積分1ステップ進める
# # if (n+1) % 100 == 0:
# # ax2.plot(x, u[nindex, 0, :], label=f"t = {dt * (n+1):.2f}")
# txt = ax.text(0, 1.01, f"t = {dt * (n + 1):.1f}", fontsize=15)
# im = ax.contourf(X, Y, u[nindex, :, :], 20, cmap=plt.cm.viridis)
# ims.append(im.collections+[txt]) # im.collectionsというところに等高線プロットの情報が入っているらしい。これと文字列を結合して、経過時間を表示している。
# cbar = fig.colorbar(im)
# cbar.set_label(r"Energy", fontsize=15)
# fig.tight_layout()
# anim = ArtistAnimation(fig, ims, interval=100)
# video = anim.to_html5_video()
# html = display.HTML(video)
# display.display(html)
# # HTML(anim.to_html5_video())
# # ax2.legend()
# plt.close(fig)
/var/folders/y5/1mt0b_g15zn1vv8_73m7hlb80000gn/T/ipykernel_79809/2864887905.py:50: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show()
u[3] = 0.9 * np.ones_like(X) + (1 - 2 * R) * (R < 0.5)
plt.contourf(X, Y, u[3])
plt.colorbar()
plt.show()
plt.plot(x, u[3,:, int(jmax-10)])
[<matplotlib.lines.Line2D at 0x13beb0bb0>]
$(x, y) -> (r, \theta)$
$$ x = r\cos \theta \\ y = r\sin \theta $$$$ e = 0.9 (r > 0.5)\\ e = 0.9 + 1 - 2*x (r < 0.5) $$はエネルギー密度